#!/usr/bin/perl
package parse_bl2seq;

require Exporter;


our @ISA = qw(Exporter);
our $VERSION = 1.00;
our @EXPORT = qw(	
					expand_cdna_with_alignment_str
					expand_cdna_with_alignment_array
					expand_prot_with_alignment_str
					expand_prot_with_alignment_array
					pos_in_aligned_seq
					get_alignment_marks
					gap_positions_array_to_str
					gap_positions_str_to_array
					);

use strict;
use warnings;

# function declarations
sub parse_blast_results(\@$$);



#_________________________________________________________________________________________

#	parse_blast_results

#		input:

#			@results
#				array of lines with the blast results

#		returns:

#			$score				bitscore
#			$alignment_length	length of the alignment
#			$i_pos				reference to an array of position / length pairs
#			$j_pos				reference to an array of position / length pairs

#_________________________________________________________________________________________
sub throw_bad_data($\@)
{
   	die "Error!\nUnexpected data on line ".$_[0]." of the blast output\n"
			.$_[1][$_[0]]."\n";
}

sub gap_positions_array_to_str(\@)
{
	return join ("|", map{$_->[0].':'.$_->[1]} @{$_[0]});
}
sub gap_positions_str_to_array($)
{
	return map{[split (/:/,$_)]} split (/\|/, $_[0]);
}
sub parse_blast_results(\@$$)
{
	my ($results, $len_a, $len_b) = @_;
	my $score;
	my ($e_value);
	my ($alignment_length, $alignment_start);
	my (@i_pos, @j_pos);
	for (my $i = 0; $i < @$results; ++$i)
	{
		# get bit score		Score =  703 bits (1814), Expect = 0.0
		next unless ($results->[$i] =~ /^\s*Score\s*=\s*(\d+(?:\.\d+)?).*Expect\s*=\s*([0-9\-\.Ee]+)/);
		$score = $1;
		$e_value = $2;
		$e_value =~ s/^e/1e/;
		
		++$i;
		++$i while (defined $results->[$i] && !length($results->[$i]));
		die "Error!\nPremature end after line $i\n" unless defined $results->[$i];

	
		# get the length of the alignment i.e. "Identities = 527/625 (84%)"
			throw_bad_data($i, @$results) unless
		$results->[$i] =~ /^\s*Identities\s*=\s*\d+\/(\d+)/;
		$alignment_length = $1;
		
		++$i;
		++$i while (defined $results->[$i] && !length($results->[$i]));
		die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
	
		# go through query/align/subject line triplets of the top alignment
		my $do_initial_align = 1;
		while ($results->[$i] =~ /^\s*Query:/)
		{
			# query
				throw_bad_data($i, @$results) unless
			$results->[$i] =~ /^\s*					# whitespace
								Query:				# name of line
								\s*                 # whitespace
								(\d+)               # start position of query
								\s*                 # whitespace
								([A-Z\-+]+)/x;		# text string
	
			my $i_aa_line	= $2;
			my $i_seq_start = $1;
			++$i;
			++$i while (defined $results->[$i] && !length($results->[$i]));
			die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
	

			#skip past alignment line
			++$i;
			++$i while (defined $results->[$i] && !length($results->[$i]));
			die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
	
	
			# subject
			throw_bad_data($i, @$results) unless
				$results->[$i] =~ /^\s*					# whitespace
									Sbjct:				# name of line
									\s*                 # whitespace
									(\d+)               # start position of query
									\s*                 # whitespace
									([A-Z\-+]+)/x;		# text string
			my $j_aa_line = $2;
			my $j_seq_start = $1;
			++$i;
			++$i while (defined $results->[$i] && !length($results->[$i]));
			die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
	
			# note the offset to the first aligned position
			if ($do_initial_align)
			{
	
				$alignment_start = $i_seq_start;
				if ($i_seq_start > $j_seq_start)
				{
					push(@j_pos, [0, $i_seq_start-$j_seq_start]);
				}
				elsif ($i_seq_start < $j_seq_start)
				{
					push(@i_pos, [0, $j_seq_start-$i_seq_start]);
					$alignment_start = $j_seq_start;
				}
				$do_initial_align = 0;
			}
	

			# save the number and position of the aligning characters for each sequence
			while ($i_aa_line =~ /(\-+)/g)
			{
				# save position and length of aligning markers
				my $len = length($1);
				my $pos = pos($i_aa_line) - $len + $i_seq_start - 1;
				if (@i_pos && $i_pos[-1][0] == $pos)
				{
					$i_pos[-1][1] += $len;
				}
				else
				{
					push(@i_pos, [$pos, $len]);
				}
				# adjust starting position to account for the markers already counted
				$i_seq_start -= $len;
			}
	
			while ($j_aa_line =~ /(\-+)/g)
			{
				# save position and length of aligning markers
				my $len = length($1);
				my $pos = pos($j_aa_line) - $len + $j_seq_start - 1;
				if (@j_pos && $j_pos[-1][0] == $pos)
				{
					$j_pos[-1][1] += $len;
				}
				else
				{
					push(@j_pos, [$pos, $len]);
				}
				# adjust starting position to account for the markers already counted
				$j_seq_start -= $len;
			}
		}
		
		if (defined $score)
		{
		#find length of two sequences
			while ($results->[$i] !~ /^length of query:\s*(\d+)/)
			{
				++$i;
				++$i while (defined $results->[$i] && !length($results->[$i]));
				die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
			}
			$len_a = $1;
			while ($results->[$i] !~ /^length of database:\s*(\d+)/)
			{
				++$i;
				++$i while (defined $results->[$i] && !length($results->[$i]));
				die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
			}
			$len_b = $1;
		}
		my $final_len_a = $len_a;
		$final_len_a += $_->[1] for (@i_pos);
		my $final_len_b = $len_b;
		$final_len_b += $_->[1] for (@j_pos);

		# add end alignment marks
		my $final_len = $final_len_a;
		if ($final_len_a > $final_len_b)
		{
			push(@j_pos, [$len_b, $final_len_a - $final_len_b]);
			$final_len = $final_len_a;
		}
		elsif ($final_len_b > $final_len_a)
		{
			push(@i_pos, [$len_a, $final_len_b - $final_len_a]);
			$final_len = $final_len_b;
		}
		return ($score, $alignment_start, $alignment_length, $e_value,
					gap_positions_array_to_str(@i_pos),
					gap_positions_array_to_str(@j_pos),
					$final_len,
					$alignment_length / $final_len);
	}
	return;

}




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	print the modified sequences

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#_________________________________________________________________________________________

#	expand_cDNA_to_alignment

#		input:

#			$seq
#			$aligning character
#			@list of position/length pairs in string format 0:63|119:6|263:1

#		returns:

#			$expanded sequences
#_________________________________________________________________________________________
sub expand_cdna_with_alignment_array($$\@)
{
	my ($seq, $char, $align) = @_;
	$char = '-' unless defined $char && length $char;
	$char = substr($char, 0,  1);
	$char = $char x 3;
	
	foreach (reverse @$align)
	{
		substr($seq, $_->[0] * 3, 0) = ($char x $_->[1]);
	}
	return $seq;

}

sub expand_cdna_with_alignment_str($$$)
{
	my ($seq, $char, $align_str) = @_;
	my @align = map {[split '\:', $_]; } split '\|', $align_str;
	return expand_cdna_with_alignment_array($seq, $char, @align)
}


#_________________________________________________________________________________________

#	expand_sequence_to_alignment

#		input:

#			$seq
#			$aligning character
#			@list of position/length pairs

#		returns:

#			$expanded sequences
#_________________________________________________________________________________________
sub expand_prot_with_alignment_array($$\@)
{
	my ($seq, $char, $align) = @_;
	$char = '-' unless defined $char && length $char;
	$char = substr($char, 0,  1);
	foreach (reverse @$align)
	{
		substr($seq, $_->[0], 0) = ($char x $_->[1]);
	}
	return $seq;
}
sub expand_prot_with_alignment_str($$$)
{
	my ($seq, $char, $align_str) = @_;
#	return unless $align_str;
	my @align = map {[split '\:', $_]; } split '\|', $align_str;
	return expand_prot_with_alignment_array($seq, $char, @align)
}

#_________________________________________________________________________________________

#	pos_in_aligned_seq

#		input:

#			$position
#			@list of alignment position/length pairs

#		returns:

#			adjusted position given the aligning marks
#_________________________________________________________________________________________
sub pos_in_aligned_seq($\@)
{
	my ($pos, $align) = @_;
	my $new_pos = $pos;
	foreach (@$align)
	{
		$new_pos += $_->[1] if ($pos >= $_->[0]);
	}
	return $new_pos;

}


#_________________________________________________________________________________________

#	get_alignment_marks

#		input:

#			$ aligned protein sequence with no spaces etc.

#		returns:

#			positions of aligning marks
#_________________________________________________________________________________________
sub get_alignment_marks($)
{
	my ($seq) = @_;
	my @align_pos;
	my $cummulative = 0;
	# save the number and position of the aligning characters for each sequence
	while ($seq =~ /(\-+)/g)
	{
		# save position and length of aligning markers
		my $len = length($1);
		push(@align_pos, [pos($seq) - $len - $cummulative, $len]);
		$cummulative += $len;
	}
	return @align_pos
}










1;

